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Abstract 

Boolean networks have been proposed as potentially useful models for genetic control. An 
important aspect of these networks is the stability of their dynamics in response to small perturba- 
tions. Previous approaches to stability have assumed uncorrelated random network structure. Real 
gene networks typically have nontrivial topology significantly different from the random network 
paradigm. In order to address such situations, we present a general method for determining the 
stability of large Boolean networks of any specified network topology and predicting their steady- 
state behavior in response to small perturbations. Additionally, we generalize to the case where 
individual genes have a distribution of 'expression biases,' and we consider non-synchronous up- 
date, as well as extension of our method to non-Boolean models in which there are more than two 
possible gene states. We find that stability is governed by the maximum eigenvalue of a modified 
adjacency matrix, and we test this result by comparison with numerical simulations. We also 
discuss the possible application of our work to experimentally inferred gene networks. 
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I. INTRODUCTION 



Boolean networks have been extensively investigated as a model for genetic control of 
cells 

HQ. 

In this model, each gene is represented by a node of a network, and each 
node has one of two states: on - i.e., producing ('expressing') its target protein - or off. 
Directed links between genes indicate that one gene influences the expression of another. 
This can correspond to the expressed protein directly binding to DNA and modulating the 
transcription of a gene or to other signaling pathways that modulate DNA transcription. In 
the standard Boolean network model, the system evolves in discrete timesteps (t = 0, 1, 2, ...), 
and at each step the state of every node is simultaneously updated according to some function 
of its inputs. This function approximates the action of activators (proteins which act to 
increase expression of a given gene) or inhibitors (proteins which act to reduce expression). 
While this model might seem to be an oversimplification considering the complex kinetics 
involved in all steps of a transcription pathway, experimental evidence suggests that real 
biological systems are, in some cases, reasonably well-approximated by Boolean networks 



In 1969, S.A. Kauffman |l[ introduced a type of Boolean network known as an N — K 
network. In this model, there are N nodes each having exactly K input links, and the nodes 
from which these input links originate are chosen randomly with uniform probability. We 
refer to the number of input (output) links to (from) a node as the in-degree (out-degree) 
of that node. At any given time t, the system state can be represented as an iV-vector 
whose ith component <Tj(i) is either zero or one, where i — 1, 2, N. There are 2 N possible 
states. The function determining the time evolution at each node is defined by a random, 
time-independent, 2^-entry truth table. Since this is a finite, deterministic system, there 
is always an attractor: eventually, the system must return to a previously visited state 
(finiteness), after which the subsequent dynamics will be the same as for the previous visit 
(determinism). These attractors can be fixed points or periodic orbits. Using the Hamming 
distance between two states (i.e., the number of nodes for which the <Ti(t) disagree) as the 
distance measure, the system exhibits both what is termed a 'chaotic' (or unstable) regime, 
where the distance between typical initially close states on average grows exponentially in 
time, as well as a stable regime, where the distance decreases exponentially. Between the 
two there is a 'critical' regime. (Here by 'close' we mean that the Hamming distance is small 
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compared to N.) 

As a model of genetic control, these attractors have been postulated to represent a specific 
pattern of protein expression which defines the cell's character [l| . In single-celled organisms, 
these attractors might be taken to correspond to different cell states (growing, dividing, 
starving, heat- or pH-shocked). In multi-cellular organisms, different cell types (muscle, 
nerve, liver, etc.) have different expression patterns, and, within each type, a cell could be 
in a variety of states (resting, 'activated,' dividing, etc.) that each correspond to different 
expression patterns. Boolean network approximations have been successful in predicting the 
gene expression time sequence of the segment polarity gene network in Drosophilia, a model 
for embryonic development where individual cells turn specific proteins on and off in patterns 
that guide the growth of certain organs and structures [3j. Since the protein expression 
pattern of the cell is modeled from the state of the corresponding Boolean network, the 
question of the stability of the network then becomes important: do small perturbations 
in the expression pattern, due perhaps to chemical fluctuations, die out quickly, returning 
the cell to its original state, or do they quickly grow, pushing the cell into another state? 
The purpose of this paper is to examine the stability of network dynamics in the context of 
discrete state models of gene networks. 

One motivation for the consideration of dynamical stability is its possible relevance to 
cancer. Specifically, we hypothesize that dynamical instability of a gene network might be 
a causal mechanism contributing to the occurence of some cancers. We emphasize that this 
hypothesis is distinct from the previous hypothesis of 'genomic instability' as a cause of 
cancer 4j. In particular, genomic instability has been defined 36] as 'the failure to transmit 
an accurate copy of the entire genome from one cell to its two daughter cells.' In contrast, 
the instability we refer to is that of the dynamics of a given gene network, and we use 
the term 'dynamical network instability' (DNI) to distinguish this condition. We speculate 
that DNI might arise from mutations and that, once established, as cells divide, DNI could 
lead to widely varying gene expression patterns from cell to cell. We emphasize that DNI 
implies that this variation would arise even in the absence of further mutation. That is, 
similar to the concept of chaos in continuous-state dynamical systems (e.g., 5]), DNI causes 
exponential sensitivity of typical system trajectories to small changes, which we speculate 
may lead to many different outcomes in the course of cell division. Recent microdissection 
results indicate wide variations in gene expression patterns even for nearby cells within the 
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same cancerous tissue g This var.abil.ty provides a basis for understanding why cancer 
can adapt and evade treatment I7fl. 

n 

Another motivation for our study is the argument, put forward by Kauffman [2[, that 
evolution favors gene networks that are on the border between stability and instability 



llO|, 111)]. Whether or not our cancer hypothesis or Kauffman's stability-border hypothesis 
holds, the question of dynamical stability of such networks is crucial to their understanding 
and use as models. 

While previous works have addressed the question of dynamical network stability in sim- 
ple, specific types of random networks (e.g. N — K nets), in this paper we address the 
question of dynamical network stability for general network topology and node attributes. 
We also consider nonsynchronous update and extend the considerations to non-Boolean 
models allowing for the possibility of nodes having more than two states. Thus our work 
provides a potentially enhanced framework for modeling and using the discrete state net- 
work paradigm. In particular, we consider how our network stability considerations can be 
employed on experimentally derived gene networks. 

In the original N — K nets as proposed by Kauffman, the truth table output governing 
node dynamics was randomly chosen with on and off having equal probability. Subsequently, 
it was shown that if the truth table output was biased such that p denotes the probability 
of randomly a ssig ning an off output, the transition between the stable and chaotic regimes 



depends on p 12]. We term p the 'expression bias.' Additionally, networks with a distribu- 



tion of in-degrees, but no in-/out-degree correlation, have been considered in [k], Jj], Jj]. 
and it has been shown that the nodal in-degree average, (K m ), suffices to determine the 
stability. (Here (•) indicates average of a nodal quantity over all nodes.) Specifically, the 
critical average number of connections, K c , governing this transition is 

K c = l/[2p(l-p)}, (1) 

where the network is stable for (K tn ) < K c , unstable for (K tn ) > K c , and critical for (K m ) = 



K c . Aldana and Cluzel [161 considered the consequences of Eq. ([I]) in the case of networks 
with scale-free topology 17|], i.e., the probability distribution P(K m ) (or P(K out )) that a 
randomly chosen node has in-degree K m (out-degree K out ) is a power-law: P(K) oc K~ y . 
(Since every out-link for a node is an in-link for some other node, (K m ) = (K out ); thus the 
result is unchanged whether it is the in- or out-degree that has power-law scaling.) 



Recently, some authors have noted, but not numerically tested, a generalization of Eq. (pQ) 
that takes into account nodal correlations between the in-degree and out-degree characterized 
by the joint K tn — K out degree distribution function P(K m , K out ). In this case, the critical 
transition occurs at [181 ] 

(K in K out ) _ 1 

(K) ~2p(l-pY () 
We emphasize that Eq. ([T]) was derived in the annealed approximation (see later discus- 
sion) for networks with a given in- or out-degree distribution P{K) and with the compli- 
mentary links completely random, and that Eq. (j2J) uses only the additional information 
contained in the nodal in- /out- degree correlation. Furthermore, all nodes ('genes') were 
taken to have the same p value. However, gene networks, in common with real networks 
occurring across a broad range of applications, can be expected to deviate substantially from 
the above simple network model. Examples of network properties that could make previ- 
ous analyses of network stability inapplicable are assortativity [3] (the tendency for highly 
connected nodes to prefer or avoid linking to other highly connected nodes) and community 



structure 20( (the existence of highly connected, sparsely interconnected subgraphs), two 
properties that are not captured in the degree distributions. Additionally, these properties 
may have biological implications. For example, a recent paper [2l| examined gene inter- 
action networks from cancerous tissue and found significant community structure, as well 
as positive correlation between the in-degree and out-degree of nodes; additionall y, p rotein 



interaction networks have been shown to exhibit significant disassortativity [19|, [22]. Fur- 
thermore, for modeling purposes, it might be important to allow the expression bias p to 
vary from node to node (as an extreme example, so-called housekeeping genes [23] have a 
predominant tendency to be on, corresponding to low p, unlike other genes). In this paper, 
we derive and test the stability criterion for large networks with arbitrary network topology 
and heterogeneous expression biases. In particular, our theory evaluates the stability of any 
given network with its specific topology (i.e., its adjacency matrix A defined subsequently), 
and by its node-specific expression biases. We show that stability is determined by the 
largest eigenvalue of a modified adjacency matrix, and we numerically test this criterion. 

With respect to real gene networks, the synchronous update at integer times (t = 
0, 1, 2, ...) used in the above models represents an additional deviation from the real situation, 
where chemical kinetics and transport processes can be expected to introduce non-trivial 
dynamics. As a partial step toward remedying this (and to make Boolean approximations 



suitable for atmospheric and geophysical processes), Ghil and Mullhaupt 24] consider a 
generalization in which t is a continuous variable and (7j(t) depends on <jj(t — Tij), where 
Tij is a delay time that can be different for each link from j to i. The original formulation 
(e.g., in Refs. 

HQ HQ) 

corresponds to r^j = 1 for all We will argue and 
numerically confirm that the criterion determining the stability/instability border of this 
generalization of the Boolean network model is the same as that for the synchronous update 
models. 

In addition to nonsynchronous update, another generalization of Boolean networks that 
we will examine is models in which each node i is allowed to have one of Si possible discrete 
states (e.g., for = 3, we label the states <jj G {0, 1,2}, and for Boolean networks Si = 2 
for all i). This model may be closer to the behavior of actual cells, and models with multiple 



states can be related to certain piece- wise ODE models of transcription |25l. |26I|. The general 
model using arbitrary, multivalued truth tables has been previously treated in the special 
case of iV — K networks with all nodes having the same number of possible states S. In the 



case where each possible state is equally likely, the critical number of inputs is [37] 



K ° = Turfs' (3) 

The applicability of our work to any specific network and set of node-wise expression 
biases may be of particular interest in situations where experimental data provide the possi- 
bility of estimating a gene network and expression biases. Such information could be used as 
input to our method which could give an indication of the stability of a given experimentally- 
derived network. The possibility that such analyses may be feasible becomes more and more 
likely with the rapid technological advances in obtaining new types of high-quality, quan- 
titative data useful for deducing gene networks. For example, such analyses could be used 
to address the hypothesis that dynamical instability of gene networks is connected with the 
occurence of cancer. 

II. MODEL 

Deterministic Boolean networks are formally defined by a state vector S(t) = 
[<7i(t)<T 2 (t)...<TAr(t)] T , where <Tj G {0, 1}, and a set of update functions fi, such that 

Vi(t) = fi(<Tk(i,l)(t - 1), (7^2) (t- I),--), (4) 



where k(i, 1), k(i, 2), k(i, Kf 1 ) denote the indices of the Kf 1 nodes that input to node i; we 
denote this set of nodes by /Q = {k(i,j)\j = 1, 2, .., Kf 1 }. The update function fa is defined 
at each node i by specifying a truth table whose outputs are randomly populated. Previous 
analytic results assumed a constant expression bias for all nodes; however, we allow that, in 
the truth table for node i, output entries are randomly assigned zero with probability pi or 
one with probability 1 — p^. In the case of uniform expression bias, we drop the subscript 
and use the notation p = pi. 

We consider the interaction structure of this system as a graph where the nodes represent 
individual elements of the state vector, and a directed edge is drawn from node j to node % 
if j G /Q. An adjacency matrix A is defined in the usual way: a matrix entry is one if 
there is a directed edge from node j to node i and zero otherwise. 

The stability of a large Boolean network is defined by considering the trajectories resulting 
from two close initial states, £(i) and To quantify their divergence, the Hamming 

distance of coding theory is used: h(t) = J2f=i \ a i{t) ~ <?i{t)\- If the network is stable, on 
average h(t) —>■ as t —>■ oo. In unstable networks, h(t) quickly increases to O(N), while a 
'critical' network is at the border separating stability and chaos. 

In order to study the stability of an N — K Boolean network, Derrida and Pomeau 



12l | considered an annealed situation and calculated the probability that, after t steps, a 
node state is the same on two trajectories that originated from init ially close conditions. 
(This calculation was later generalized to variable in-degree [l3, 14, 15j, and joint degree 
distribution in [la].) In Derrida and Pomeau's 'annealed' situation, at each time step t 
the truth table outputs and the network of connections are randomly chosen. The actual 
situation of interest, however, is the case of 'frozen-in' networks, where the truth table and 
network of connections are fixed in time. It has been commonly assumed that analytical 
results obtained for the annealed case are a good approximation to the frozen-in case (e.g., 
Refs. [12I, 13, Q, [15I). We also adopt this view in a modified form, and we will test its 
predictions with numerical simulations. 

The randomization of the network of connections at each time step while keeping the 
degree distribution fixed carries the implicit assumption that there is no additional dy- 
namically relevant structure in the frozen network other than that contained in the joint 
degree-distribution P(K m , K out ). To avoid this assumption, we obtain theoretical results 
for a different annealing protocol, which we term 'semi-annealed.' In this semi-annealed pro- 
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cedure, we keep the network fixed (i.e., the adjacency matrix A does not change with time), 
and we envision randomly assigning the output entries of the truth table of each node i at ev- 
ery time t according to the time-independent expression bias Pi assigned to node i. We then 
imagine tracking the probability that individual node states cij(t) and (Xj(t) differ over time 
with an iV-dimensional difference vector, whose components are Ui{t) = ((\<Ji(t) — <fj(i)|)), 
where ((•)) denotes an average over every possible small initial perturbation. Here by 'every 
possible small initial perturbation' we mean all perturbations for which a small fraction e 
of the states are flipped. Additionally, we define the 'sensitivity' qi as the probability that 
the output of fi changes when given two different input strings, similar to the 'average 



sensitivity' of Ref . [27| . In the case of completely random Boolean functions, 

ft = 1 - (Pi + (1 - Pi) 2 ] = 2p<(l - Pi)- (5) 



Thus, similar to Ref. [121 ] . we can write the update equation for i/i as 



y t (t) = q l (l-H(l-y J (t-l))\ (6) 

^ je/Q ' 

Equation ([H]) follows from noting that the probability that <Jj and dj are equal is (1 — yj) 
and thus the probability that all inputs to node i are equal is the above product. Note that 
this equation uses topological information contained in the /Q. However, we have treated 
Uj, yji and qi as if they were probabilities of statistically independent random events. We 
hypothesize that this semi-annealed protocol might be expected to yield good results for 
frozen-in cases when the network is large and the fraction of network nodes on short loops is 
small (the network is 'locally tree-like'). To see the problem posed by short loops, consider a 
node with two inputs that themselves have inputs both coming from a common node; in this 
case, the elements of y(t) in Eq. dSJ) are no lon ger statistically independent and multiplying 



28| for discussion related to the locally tree- 



the probabilities is no longer correct. See Ref. 
like assumption. Our numerical tests of frozen networks indeed yield results that agree very 
well with our semi-annealed hypothesis on large, locally tree-like networks. We also find our 
predictions to hold for networks with a large number of feedforward motifs, a nontree-like 
three-node subgraph that has been found to be prevalent in real gene networks 30]. 



The case where both network states are exactly the same corresponds to Ui{t) = 0, which 
is a fixed point of Eq. ([6]). In order to determine the stability of this fixed point, we linearize 



S 



Eq. dSJ) around y(t) = for small perturbations: 

N 

Vi(t + 1) « gi^A^, (7) 

where Ay are the elements of the adjacency matrix A. Equation (j7j) can be written in matrix 
form as y(t + 1) = Qy(t) where 

Qij = qiAij. (8) 
The stability is thus governed by the largest eigenvalue Xq of this matrix: 

\q > 1, y = is unstable; 

Xq = 1, y = is critical; (9) 
\q < 1, y = is stable. 



Since Qy > 0, the Perron-Frobenius theorem [3JJ guarantees that Aq is real and positive. 
We also note that, for any given adjacency matrix A and assignment of q^s to nodes, Eq. 
(EJ) can be iterated numerically to predict the expected time-asymptotic saturation value of 
the difference in two initially nearby states when evolved to steady-state. We numerically 
test this prediction, as well as the stability criterion in Eq. ([9]) in the next section. [3^] 

As a special case of interest, if the the g« are uniform, qi = q, then Xq = qX, where A is 
the maximum eigenvalue of the adjacency matrix. This yields the critical condition, 

A = 1/q. (10) 

Furthermore, for the case of a large network whose links are randomly assigned subject to 
a joint probability distribution P(K tn , K out ) at each node (with no assort at ivity), the mean 
field approximation for the largest eigenvalue is 3^] 

/ l/'in Tsout\ 

X pa '- (11) 

where, since (K m ) = (K out ) necessarily, we use the notation (K) = (K m ) = (K out ). Equa- 
tions (TTOT) and ( TlTl) yield the same criterion as in Eq. ([2j). In the case where K m and K out 
are uncorrelated, P(K in ,K out ) = P m (K m )P out (K out ) and (K in K out ) = (K) 2 , yielding Eq. 

©■ 

The eigenvalue of random network adjacency matrices with assortativity has been con- 



sidered in Ref. 



321 ]. which defines an assortativity measure p as 



(Kl n K° ut ) e t s 

n=- - 3 - — (12) 

(K in K out ) ' 



where (K\ n K° ut ) e denotes an average over all links (i, j) from node i to node j. The network 
is assortative (disassortative) if p > 1 (p < 1). For p near one, the largest eigenvalue A is 



approximately given by [32] 



/ jsin jy~out\ 

^ ~ (K) ' < 13 » 
Thus by Eqs. ([6]) and ([9]), it is predicted that, for uniform q, assortativity (disassortativity) 
decreases (increases) the critical q value. 

In the case of nonuniform g i; we have recently generalized Eq. ffTTT) to obtain an analogous 
mean field approximation to \q without assortativity or community structure, 

( q K*"K°«) 

Aq ^— - (14) 

Our derivation of (jbil) will be published elsewhere. From ffT4l) . we see that correlation (anti- 
correlation) between q and K m K out decreases (increases) network stability and that, in the 
absence of correlation, the result is similar to that for a uniform q, Xq « (q) (K %n K owt ) / (K), 
with (q) replacing the uniform q (Eqs. ([9]) and (|T0l)). 

We now consider the generalization to allow any number of discrete node states. We 
denote the number of possible states of node i by Si, and we label the possible states 
0, 1, 2, Si — 1. The number of possible inputs to i from the set /Q of nodes that influence 
it is rijeyCi Sj- For each of these possible inputs, the truth-table function f { in Eq. (jlj) assigns 
one of the Si possible states to node i. Similar to the Boolean case, we take the assignment 
to be random and to have an 'expression bias' pi^ s for each of the s — 0, 1, 2, ...Si — 1 node 
states, where p i)S denotes the probability that f iy for a given set of inputs, assigns the state 
s to node i, and Yl s Pi,s = 1 f° r an n °des i. As in the Boolean case, we can then introduce 
the sensitivity q,i giving the probability that two different sets of inputs result in a different 
updated state of node i, which, in the random truth table case, 

ft = 1 (15) 

s 

With this definition, we see that all our previous reasoning still applies, and Eqs. ([6])- (Q 
hold with this generalized expression for the node sensitivities and with Di{t) interpreted as 
the probability of disagreement between <7j(t) and di(t). In the case of uniform number of 
node states Si = S and equal expression biases pi )S = p s = 1/S, among these states, Eq. 
( TL5]) becomes q = 1 — 1/S, which, when combined with Eq. (fTUl) yields the previous result 
inEq. ©. 
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Finally, we note that our criticality criterion, Xq = 1, is unchanged by the presence of 



delays, as in the models of Refs. 24| , and only a slight modification is required of Eq. ([6]) 
(i.e., Ujit — Tij) replaces Uj(t — 1)). The condition Xq = 1 implies that the components of y in 
Eq. (jnD are time-independent. Thus we predict that the delays Ty do not influence the result, 
and the criticality condition in Eq. is independent of the synchronous update structure of 
the most commonly used random Boolean network models. Similarly, the time-asymptotic 
steady state obtained by repeated iteration of §6§ is, by definition, time-independent and 
thus also does not depend on the r^- (although the will influence the time-dependent 
approach to the asymptotic steady state; see supplementary material). 

III. STATISTICAL METHODS 

We numerically test the above predictions on several classes of Boolean networks with 
uniform sensitivity (i.e., = q is the same for all nodes): 

(a) random networks with K m = K out ; 

(b) random networks with imperfect correlation between K m and K out ; 

(c) networks with assortativity or disassortativity; and 

(d) networks constructed as in (a) but with a substantial number of feedforward loops. 

We additionally test our predictions on two classes of networks with nonuniform sensi- 
tivities: 

(e) networks constructed as in (a) but where nodes have different sensitivities correlated 
with the degrees of the nodes; and 

(f) networks with significant community structure, where the two communities have dif- 
ferent, uniform sensitivities. 

Finally, we test our generalization to more than two node states on networks of type (a) 
but with Si = 4 for all nodes. For types (a)-(c) and (e), we use networks with truncated 
power-law degree distributions. (Evidence for the presence of this type of distribution in 



gene networks has been seen in 33].) 
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The algorithms for constructing the networks of types (a)-(c) are as follows, (i) Establish 
the in- and out-degrees for each node, which are drawn from a distribution, 

{K~ 7 , K < K max } 
(16) 
0, K > K max , 

where 7 = 2.1 and K max = 15 (Boolean case) or K max = 8 (Si = 4 case). The out- 
degree is initially set to the in-degree. (ii) Randomly swap the out-degrees between pairs 
of nodes. If maximal correlation between in- and out-degrees is desired, as in (a), this 
step is skipped so that K m = K out and (X in K out ) is maximal. A completely uncorrelated 
network has every nodal out-degree swapped exactly once, yielding ^K in K out ) = (K) 2 . The 
quantity (K m K out ), which approximately determines A by Eq. ffTTl) . can thus be tuned by 
the number of nodes that have their out-degrees swapped, (iii) Place links randomly between 
nodes subject to the constraints of the specified in- and out-degrees assigned at each node 
by the 'configuration model' 



34j . (iv) If networks with assortativity (disassortativity) are 



desired, as in (c), perform a given number of link swaps, as in [32j, that increase (decrease) 
the assortativity p in Eq. ffl2|) . In all cases we employ networks with iV = 10 4 and two 
initial conditions separated by a Hamming distance of 100. In the supplementary material 
we discuss finite size effects that can occur for smaller N. 

We emphasize that, although we determine our networks randomly, in our numerical 
experiments we do not average over this randomness. Rather, we generate one random 
network for each experiment and examine the resulting behavior of that specific network. 



IV. RESULTS 



We test the steady-state predictions of Eq. and the criticality condition of Eq. ([9]) in 
Fig. [H and compare the calculated critical parameters to the mean-field-type approximations 
of Eqs. p3D , (fT2"j) . and ffT4l) in the supplementary notes. In order to compare Eq. §6§ (solid 
curves in Fig. [1]) to experimental measurements of the Hamming distance from numerical 
evolution of true frozen Boolean dynamical systems (markers in Fig. [1]), we calculate the 
node averaged steady-state fractional Hamming distance, 

y=&^Ev*(*)- ( 17 ) 
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In practice, this limit is calculated as the average Hamming distance from t = 90 to t = 100 
when all delays are the same (r^ = 1), and from t = 490 to t — 500 when nonuniform delays 
are present. These times are well after the steady-state value is reached (see supplementary 
material). Each experimental data point in Fig. [TJ corresponds to a single realization of 
interconnections averaged over 100 realizations of the time-independent truth table with 
specified sensitivity as before. 

A. In- /Out- degree Correlations and Heterogeneous Time Delay 

Figure QJa) shows the steady-state Hamming distance as a function of the sensitivity for 
one network of type (a) (A = 4.4) and two of type (b) (A = 2.9,2.3). The closed markers 
in the figure represent experiments with uniform delay Ty = 1 on all links, while the open 
markers correspond to experiments where half the links, randomly chosen, have Ty = 10 
and the remainder have Ty = 1. Importantly, the degree distributions are the same for 
all three networks, and we attain different A values by varying the correlation between the 
in-degree and the out-degrees. We see from Fig. Ufa) that there is close agreement between 
the theoretical prediction and the experimental results and our prediction that the presence 
of delays does not change the stability is confirmed. Additionally, the measured steady-state 
Hamming distance is essentially zero below the critical value of the sensitivity, q cr u = 1/A 
(this point is indicated by vertical downward arrows in Fig. [TJ^a)). We emphasize that the 
degree distributions (and hence (A")) are the same for the networks in Fig. [H(a), and thus, 
if the in-/out-degree correlation were ignored, the observed difference between the stability 
conditions for these networks would not be predicted. 

B. Assortativity/Disassortativity 

Figure d^b) shows results obtained when significant assortativity or disassortativity is 
present (type (c) networks). In this experiment, as well as all those reported below, the delays 
are all uniform. The networks under consideration have the same joint degree-distribution 
with K m = K out . However, each of the networks have very different assortativities (p = 
0.52, 1.0, 1.7, defined in Eq. (fT2]) ). which yield different largest eigenvalues (A = 3.0, 4.4, 9.9). 
Since the joint degree distributions are the same, Eq. (j2J) would predict that the three 
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networks have the same stability characteristics. However, since their eigenvalues are very 
different, we predict that, as observed, the transitions of the three networks occur at different 
values of q. Again the theoretical predictions of q crit are indicated by vertical arrows. 

C. Motifs 

Random construction of networks, as used in the networks above, is expected to yield 
networks that are locally tree- like [32|]. However, we note that biological and other types 
of networks often have motifs (small subgraphs) that occur with higher frequency than in 
randomly constructed networks [3^]. For gene networks of E. coli and S. cerevisiae, it was 
found that the number of feedforward loop motifs (see inset to Fig. [H(c)) is significantly 
enhanced compared to the expected number in a randomly constructed network. In these 
networks, the number of feedforward loops per node c is roughly 0.1. Thus we consider a 
network of type (a) (A ~ 2.9 and N = 10 4 ) after adding 1000 (c = 0.1) and, in an extreme 
case, 2000 (c = 0.2) feedfoward loops. To add a feedforward loop, we randomly choose 
a node A, follow a random output to node B, and follow a random output of B to node 
C . We then add a link from node A to node C. We do this a given number of times, 
avoiding nodes that already participate in an added feedforward loop. In Fig. [U(c), we see 
that the semi-annealed theory of Eq. (J6l) (solid curve) again agrees well with our numerical 
experiments (solid markers). Based on such results, we believe that the locally tree-like 
network requirement does not invalidate application of our method to real gene networks. 
We also note that the critical point is essentially unchanged by the addition of loops (adding 
links only slightly increases the largest eigenvalue), however more feedforward loops tend to 
increase the steady-state Hamming distance for q > q cr u- 

D. Application to S. cerevisiae 

As a real biological example, we include in the supplementary notes a graph similar to 
those in Figs. [11(a)- (c) using a published network for the yeast S. cerevisiae [35(. 
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E. Heterogeneous Correlated Sensitivities 

Figure QJd) demonstrates the effect of a distribution of q^s on the stability of a network 
with K\ n = K° ut = Ki and with correlation between the nodal values of g^ and i.e., type 
(d) networks. We consider two situations, one where (qK 2 ) / ((q) (K 2 )) is maximal, and one 
where it is minimal. The qi are drawn from a uniform distribution centered at go (the abscissa 
in the figure), with width Aq = 0.1. Maximal (minimal) (qK 2 ) is attained by assigning the 
largest qi to the node with the largest (smallest) K i: the second largest qi to the node with 
the second largest (second smallest) Ki, and so on. As can be seen from the figure, there is 
good agreement between the semi-annealed theory and the numerical experiments, and the 
two networks become unstable at different values of go- (Vertical arrows again indicate the 
points where Aq = 1.) 

F. Community Structure 

Figure [T^e) shows our results for a case where there is community structure and 
community-dependent sensitivity. To construct the networks in Fig. [H(e), consider the 
case where there are two communities, and we assign a link from node i in community a 
to node j in community b with probability 9^. We impose the additional constraints that 
daa = &bb = @u and that 9 a b = &ba = On, and the size of the two communities are the same, 
N/2. We take (K in ) = (K out ) = (K) = (6 U + 6 n )N to be the same for both communities, 
and we also assume that communities a and b have different sensitivities q a and g&, respec- 
tively. As 8 n is increased from zero to 6 n = U; Aq changes from the case of two completely 
separated communities to one of a single random network. Communities a and b have equal 
sizes of 5000 nodes, community a has q a = 0.5, and community b has g& = 0.1. In order to 
vary Aq, we vary 8 U and 8 n , keeping their sum constant in order to maintain constant (K). 
As with the curves in Fig. [T](a)-(c), the transition to chaos is governed by Aq (Aq = 1 at 
the vertical arrow), and Eq. (J6]) (solid curve) accurately predicts the numerically observed 
(solid circles) steady-state Hamming distance. 
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G. Non-Boolean Models 

Figure [TJf) illustrates an application to a case in which there are more than two possible 
states at each node. In particular, we consider S = Si — 4 possible states at each node. 
(Since the number of possible inputs to the truth table for node i in this case is 4 I", we take 
K max = 8 due to memory constraints.) Labeling the possible node states Oi G {0, 1,2,3}, 
we take nodes to have uniform expression biases for occurence of state-label 0, po = p^o, 
from to 1. The three remaining labels (a = 1, 2, 3) also have uniform biases for all nodes 
h Ps = Pi,s = (1 - Pi,o)/3. From Eq. ([15]), q = q { = 1 - \p\ + (1 - p Q ) 2 /3], which has a 
maximum q max = 0.75 at p = 0.25. As can be seen in the figure, the predicted fraction of 
nodes with differing states y (solid curve) also has a maximum there. It is again seen that 
the measurements (markers) are well-predicted by the theory. 

V. DISCUSSION 

In this paper, we have presented theoretical results (Eqs. ([6]) and Q) which predict the 
steady-state Hamming distance between states evolved from two nearby initial conditions 
and the stability of a given network. These results are derived using the hypothesis that a 
theory derived in the semi-annealed case approximates the true situation, where by semi- 
annealed we mean that the network of connections is frozen, but the truth table at each node 
is randomly reassigned at each timestep. For large networks, this approximation was found 
to give excellent agreement with the true case of frozen connections and frozen truth tables. 
Our semi-annealed hypothesis does not rely on gross statistical properties of the network, 
but instead uses the specific network topology, as characterized by the network adjacency 
matrix, and the individual node sensitivities to make predictions. 

We tested our theoretical predictions with numerical experiments. Previously unad- 
dressed issues that we considered include the effects of assort at ivity, nonuniform time delay, 
nonuniform sensitivity, motifs, and community structure. In all cases tested we found good 
agreement with our theory. 

The theory that we have presented and tested above may represent a step forward in 
facilitating the application of discrete state dynamical network models to biological systems. 
Given a specific genetic interaction network and an estimate of the node sensitivities, Eq. ([9]) 
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FIG. 1: (a) y vs. q for three networks with different largest eigenvalues (A ~ 5.5,3.4,2.3), both 
with uniform delay on all links Tjj = 1 (closed markers) and with half the links having increased 
delay of Tjj = 10 (open markers). The solid curves correspond to the prediction y (defined in Eq. 
(|17p ) obtained by simulating Eq. ©. The downward vertical arrows correspond to q cr u = 1/A for 
each of the three networks, (b) y vs. q for three networks with different assortativities. (c) y vs. q 
for networks with added feedforward motifs, with an illustration of the feedforward motif (inset), 
(d) y vs. qo for networks with maximum correlation (circles) and minimum correlation (squares) 
between K\ n K° ut and qt, where q% is drawn from a uniform distribution centered at qo with width 
0.1. (e) y vs. #n/(#u + $n) for networks with community structure, where the two communities 
have q a = 0.5 and qb = 0.1. (f) y vs. po for a network where each node can take one of Si = 4 
possible states, pq is the probability that a zero appears in the truth table output; the remaining 
three symbols appear with equal probability. 

predicts the stability of that particular network directly from the adjacency matrix. Curated 
networks already exist in the literature for model single-cellular systems, and new algorithms 
continue to be developed for inferring interaction networks from a wide range of data sources 
(microarray experiments, GO annotation, genome sequencing, etc.). We note that such a 
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procedure has the advantage that, because the actual experimentally determined network 
is employed, topological aspects such as nodal in-/out- degree correlation, assort at ivity, 
community structure, etc., do not first have to be determined and then statistically modeled. 
Thus, by use of our stability criterion ([9]), there is the potential that future analysis may 
be able to evaluate a supposed relationship between the stability characteristics of various 
networks and their functioning. For example, one might test whether cancer gene networks 
are less stable than those in healthy tissue. This could lead to the strong variations in 
gene expression observed in cancerous tissue jf|, even when the underlying gene network is 
unchanged. We are currently pursuing research along this line. 

This work was supported by NSF (Physics) and by ONR (contract N00014-07- 1-0734). 
We thank L. Staudt for discussion. The work of A. P. was partly supported by the NCI 
intramural program. 

VI. SUPPLEMENTARY INFORMATION 1: TRANSIENT EVOLUTION 

Figure [Si] shows the time evolution results for networks of type (a) (i.e., Kf 1 = K° ut for 
all i and uniform q% = q). Once a random network is generated, we simulate the evolution 
of two close initial conditions and plot the Hamming distance as a function of time in Fig. 
ISTl Specifically, we take an arbitrary initial condition and generate a perturbed initial 
condition by flipping a fraction e of the state bits; in Fig. ISTl e = 0.01 and N = 1000 
corresponding to 10 flipped bits. Figure ISTTa) shows the Hamming distance as a function 
of time step t for four cases with different values of sensitivity (q = 0.5, 0.4, 0.3, 0.2) and 
uniform delays r%j = 1 (as in Eq. [4]). Each of these four curves are generated using the 
same network of interconnections (for which A = 4.3) and the same perturbation in initial 
conditions averaged over 100 realizations of the nodal truth tables. For the three cases 
q = 0.5, 0.4, 0.3, \q = q\ > 1, the network is predicted to be unstable. We see in Fig. ISlT a) 
that in these cases, the Hamming distance rises and eventually saturates at a constant value. 
In the fourth case, q = 0.215 and we have that qX < 1, and the network is predicted to 
be stable, which is demonstrated in the figure. These cases illustrate the strong effect that 
in-/out-degree correlations can have: for the value (K) = 1.89 in the network of Fig. ISTl 
the prediction from the result for uncorrelated networks, Eq. [1], is stability for all values 
of q (the minimum value of 1/q, the right hand side of [1], is 2, which exceeds (K) = 1.89). 
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FIG. SI: (a) Evolution of the Hamming distance between two initial conditions for a typical network 
of size N = 1000 and e = 0.01 for various values of the sensitivity and uniform delay Tij = 1. (b) 
Evolution of the Hamming distance between two initial conditions for a typical network of size 
N = 1000, q = 0.5, and e = 0.01. Results are shown for a network with Tj,- = 1 for all links (solid 
curve) and with T%j = 10 on 0.1 (dashed curve) and 0.5 (dotted curve) of the links. 

Figure ISl( b) shows time traces of the Hamming distance for q = 0.5 when non- uniform 
delays are present. In the curves shown, a fraction T = 0, 0.1, 0.5 of the links are randomly 
chosen and given delays of = 10 with the remaining links having delay = 1. The 
curves are for the same network as in Fig. [STT a) with q = 0.5 and are again the average 
of 100 different realization of the truth table. (The choice of delayed links is the same for 
all 100 realizations.) In each case, we see that the network is unstable and the Hamming 
distance rises to the same steady-state value, albeit at a slower rate for larger T. This result 
thus is consistent with our prediction that whether or not a network is stable and its final 
saturation value do not depend on heterogeneity of the delays. 

VII. SUPPLEMENTARY INFORMATION 2: FINITE-SIZE EFFECTS 

In Fig. IS2T a) and (b), we consider the importance of finite-size effects by varying e (a 
parameter which does not appear in the theory) for two different size networks of type 
(a). Figure IS2l(a) also compares the results of simulating the frozen case (solid markers) to 
the semi-annealed case described in the Theory section (open markers) for N = 10 3 . As 
before, in simulating a semi-annealed network, at each time step the nodal truth tables are 
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randomly generated with the same q. The networks under consideration in Figs. 182( a)- (b) 
have T = 0, A ~ 5.0, and (a) N = 10 3 and (b) iV = 10 4 . As the figure demonstrates, larger e 
yields better agreement with the theory, and the semi-annealed case seems indistinguishable 
from the frozen case. Note also that the results for e = 0.01 (e = 10 -3 ) and N = 10 3 is 
similar to that for e = 10~ 3 (e = 10~ 4 ) and iV = 10 4 , suggesting that the relevant quantity 
is eN, the number of flipped states. The inset of Fig. 1521(a) shows the histogram of the 
Hamming distances used in calculating the point q = 0.4, e = 0.01 (upward vertical arrow). 
The different trials used in generating this histogram correspond to different truth table 
realizations. The distribution shown in the inset consists of a large number of samples with 
Hamming distance zero and a roughly symmetric part that has a mean near the theoretical 
prediction. The overestimation of the mean by the theory therefore seems to be driven by 
the relative number of zero samples compared to the symmetric part. 

In order to understand the origin of the zero samples, we note that one way that they 
can arise is through 'irrelevant' nodes (i.e., nodes that do not influence the dynamics of the 
network) and 'frozen' nodes (i.e., nodes whose output is independent of its inputs due to the 
random assignment of the truth table). Irrelevant nodes can arise by either having no out- 
going links or by inputting only to other irrelevant or frozen nodes. Flipping the value of an 
irrelevant node, by definition, does not change the subsequent evolution of the network; if a 
perturbation between nearby initial conditions consists solely of such flips, that perturbation 
dies out quickly. Assuming that the fraction of irrelevant nodes is independent of N, then 
the probability that all eN nodes for which the two initial conditions differ are irrelevant 
goes to zero as N — > oo for constant e; in this case, the observations should exactly match 
the theoretical prediction. This is consistent with the trend indicated by our comparison of 
the N = 10 3 network in Fig. [S2|a) with the N = 10 4 network in Fig. 152Tb) . 

VIII. SUPPLEMENTARY INFORMATION 3: COMPARISON OF MEAN-FIELD 
EIGENVALUE APPROXIMATIONS WITH EXACT \ Q 

For the systems tested in our numerical experiments, Table [51] shows the critical pa- 
rameter values at the stability/instability border as obtained from direct calculation of the 
maximum eigenvalue of the matrix Q for the relevant specific networks (downward arrows 
in Fig. 1 (a)-(f)) compared to the corresponding results predicted from the mean-field-type 
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FIG. S2: The steady-state fractional Hamming distance h/N for (a) N = 10 3 and (b) N = 10 4 as 
a function of the sensitivity q for various values of e, both in the frozen case (filled symbols) and 
the annealed case (open symbols). The largest eigenvalue of this network's adjacency matrix is 
A m 5. While the theory does not depend on the value of e, finite-size effects cause a dependence 
on the number of flipped bits. The inset to (a) shows a histogram of measured Hamming distances 
at q = 0.4 and e = 0.01 (up arrow). 

theoretical estimates (Eqs. [11], [13]- [14], and [18]). For the community structure example 
we use the approximation, 



which applies for the case of two equal communities with symmetric connectivity probabil- 
ities (9 a b = Oba = On,@aa = &bb = #u) as i n Fig- l( e )- The analysis leading to ( |T8l) will be 
published elsewhere. 

The largest eigenvalue approximations predict the observed transition to unstable behav- 
ior quite well, as seen in the table below. The only exceptions to this agreement are in the 
case of significant assortativity or disassortativity; however, this is to be expected since the 
approximate theory is a linear approximation for values of p close to one. The values of 
assortativity and disassortativity used in the paper (1.7 and 0.5) are far from this regime. 
Nevertheless, even for these cases, the theory correctly predicts the qualitative trend that 
assortativity (disassortativity) decreases (increases) the critical q. 




0u(q a + q b ) + [OuiQa -Qb) 2 + ^0 2 n q a q b ] 



(18) 
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TABLE SI: Comparison between criticality conditions evaluated directly from Q and from the 
approximate theory. 





Direct Evaluation Approximate Theory 


Critical q J s from Fig. 1(a) 


n 99 


U.ZO 


(Approx. Theory from Eq. [Ill) 




n id 

u. o^t 




U.4o 


n /i/i 

U.44 


Critical q's for Fig. 1(b) 


U. 1U 


U. lO 


(Approx. Theory from Eq. [121) 


U.zz 


U.zo 




u.oo 


D 4^ 


Critical q's for Fig. 1(c) 


U.oo 


U.oo 


(Approx. Theory from Eq. [Ill) 


fl 1A 


u.oo 




n ia 


n 34 


Critical an for Fig. 1(d) 


0.19 


0.20 


(Approx. Theory from Eq. [14]) 


0.27 


0.28 


Critical 6> n /(0u + Or) for Fig. 1(e) 


0.21 


0.21 


(Approx. Theory from Eq. [18]) 






Critical po for Fig. 1(f) 


0.80 


0.80 


(Approx. Theory from Eq. [15]) 







IX. SUPPLEMENTARY INFORMATION 4: APPLICATION TO THE REGULA- 
TORY NETWORK OF S. CEREVISIAE 

Figure IS3l similar to Figs. l(a)-(c), illustrates an application of Eq. [17] to the published 
network of the yeast S. cerevisiae [35]. The largest eigenvalue of this network is A = 2.5. We 
have assumed in this plot that each node has the same sensitivity q, and again we see that the 
network undergoes a transition from stable to unstable behavior at q cr u = 1/A. However, in 
order to draw any conclusions about the criticality of the yeast regulatory network, we need 
a reliable estimate of the individual g,'s, from which we can calculate Xq. While estimating 
the sensitivities may be possible with existing microarray datasets, this is beyond the scope 
of this paper. 
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FIG. S3: y vs. q, calculated from Eq. [17], for the published regulatory network of S. cerevisiae 
[35]. The network undergoes a transition from stable to unstable behavior at q cr it = 1/A = 0.40. 
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percolation on directed networks where different sites have different removal properties. A 
similar condition involving (K 2 )/(K) also arises in percolation on undirected networks 29]. 
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